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A level-set method is developed for numerically capturing the equilibrium solute-solvent interface 
that is defined by the recently proposed variational implicit solvent model (Dzubiella, Swanson, 
and McCammon, Phys. Rev. Lett. 104, 527 (2006) and J. Chem. Phys. 124, 084905 (2006)). In 
the level-set method, a possible solute-solvent interface is represented by the zero level-set (i.e., the 
zero level surface) of a level-set function and is eventually evolved into the equilibrium solute-solvent 
interface. The evolution law is determined by minimization of a solvation free energy functional that 
couples both the interfacial energy and the van der Waals type solute-solvent interaction energy. The 
surface evolution is thus an energy minimizing process, and the equilibrium solute-solvent interface 
is an output of this process. The method is implemented and applied to the solvation of nonpolar 
molecules such as two xenon atoms, two parallel paraffin plates, helical alkane chains, and a single 
fuUerene Ceo- The level-set solutions show good agreement for the solvation energies when compared 
to available molecular dynamics simulations. In particular, the method captures solvent dewetting 
(nanobubble formation) and quantitatively describes the interaction in the strongly hydrophobic 
plate system. 

PACS numbers: 61.20.Ja, 68.03-g., 82.20.Wt, 87.16.Ac, 87.16.Uv 



I. INTRODUCTION 

The correct description of solvation free ener- 
gies and detailed solution structures of biomolecules 
is crucial to our understanding of molecular pro- 
cesses in biological systems. Efficient theoretical ap- 
proaches to such descriptions are typically given by 
so-called implicit (or continuum) solvent models of 
the aqueous environment [H, In those models, 
the solvent molecules and ions (e.g., as in physio- 
logical electrolyte solutions) are treated implicitly 
and their effects are coarse-grained. In particular, 
the description of the solvent is reduced to that of 
the continuum solute-solvent interface and related 
macroscopic quantities, such as the surface tension 
and the position-dependent dielectric constant serv- 



*e-mail address: lcheng@math.ucsd.edu| 
t e-mail address: jdzubiel@ph.tum.de 
t e-mail address : j mccammonOucsd. edu] 
5e-mail address: bli@math.ucsd.edu 



ing as input or fitting parameters. 

Most of the existing implicit solvent models are 
built upon the concept of solvent accessible surface 
area (SASA) defined in several ways [1, In these 
models, the solvation free energy is proportional to 
the SASA for the nonpolar contribution, comple- 
mented by the Poisson-Boltzmann (PB) ^, Q or 
Generalized Born (GB) 8, 9] description of electro- 
statics, i.e., the polar contribution. Although suc- 
cessful in many cases, the general applicability of 
these rather empirical models with many system- 
dependent, adjustable parameters (e.g., individual 
atomic surface tensions) is often questionable, when 
compared to more accurate but computationally ex- 
pensive explicit molecular dynamics (MD) simula- 
tions or experimental results. It is believed that the 
key issues here are the decoupling and separate anal- 
ysis of surface area, dispersion and polar parts of the 
free energy, and the inaccurate free energy estima- 
tion due to a predefined solvent-solute interface, an 
ad hoc input. It is additionally well established by 
now that cavitation free energies do not scale with 
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surface area for high curvatures [13, [ill, a fact of 
critical importance in the imphcit modehng of hy- 
drophobic interactions in biomolecular systems [12l | . 

Recently, Dzubiella, Swanson, and McCammon 
, [l^ have developed a variational implicit solvent 
model. The basic idea of this approach is to intro- 
duce a free energy functional of all possible solute- 
solvent interfaces, coupling both the nonpolar and 
polar contributions of the system, and allowing for 
curvature correction of the surface tension to ap- 
proximate the length-scale dependence of molecular 
hydration. Minimizing the functional leads to a par- 
tial differential equation whose solution determines 
the equilibrium solute-solvent interface and the min- 
imum free energy of the solvated system. This stable 
solute-solvent interface is an output of the theory. 
It results automatically from balancing the different 
contributions of the free energy. First applications 
of this approach to simple, highly symmetrical so- 
lutes showed promising results when compared to 
MD simulations 13,, 14]. 

In this work, we develop a level-set method for nu- 
merically capturing arbitrarily shaped solute-solvent 
interfaces determined by the solvation free energy 
functional in the variational implicit solvent model. 
In our method, a possible solute-solvent interface 
is represented by the zero level-set (i.e., the zero 
level surface) of a level-set function; and an ini- 
tial surface is evolved eventually into an equilibrium 
solute-solvent interface. The level-set method is a 
general technique for numerically tracking moving 
fronts with possible topological changes such as in- 
terface merging and break-ups [H, [la, [13 • Previ- 
ous applications include two-phase fluid flow, crys- 
tal growth, materials modeling, shape optimization, 
imaging process and graphics, etc.; see [H, [13 for 
related references. Recently, Can, Chen, and Wang 
[ist used the level-set method for imaging a large 
biomolecular surface based on a SASA type model. 
Our work is quite different: we not only represent 
molecular surfaces using level-sets, but further de- 
velop an evolution algorithm that numerically deter- 
mines the stable equilibrium solute-solvent interface 
based on physical rationale. 

Our new level-set techniques include two efficient 
and stable methods. One is a careful treatment of 
singularities formed during the level-set evolution 
based on certain geometrical motion of the molec- 
ular surface. The other is a two-grid numerical 
method for the calculation of the free energy using 
a Lennard- Jones type potential which changes dra- 
matically for short range interactions. 



We apply our method to the solvation of nonpo- 
lar molecules such as two xenon atoms, two paraf- 
fin plates, model helical alkane chains, and a single 
Ceo molecule. Our extensive numerical results show 
good agreement with MD calculations. In particu- 
lar, our method is able to capture the solvent dewet- 
ting phenomenon, i.e., formation of a nanobubble 
within the strong hydrophobic confinement caused 
by the paraffin plate arrangement. Furthermore, we 
demonstrate that topological changes, such as the 
rupture of a nanobubble and the fusion/breakup of 
the surface of two molecules, are captured by our 
level-set method. 

The rest of the paper is organized as follows: In 
Section ini we review the variational implicit solvent 
model. This is followed by a description of our level- 
set method in Section Hill In Section [TVl we report 
the numerical results of our level-set calculation of 
several nonpolar systems. Finally, in Section |Vj we 
conclude and give an outlook to further necessary 
extensions of our approach. 



II. VARIATIONAL IMPLICIT SOLVENT 
APPROACH 

In the following, we briefly summarize, with a few 
remarks, the variational implicit solvent approach 
that has been recently proposed by Dzubiella, Swan- 
son, and McCammon |13l . [l4| . 



A. Geometry 

Consider the system of an assembly of solutes with 
arbitrary shape and composition surrounded by a 
solvent. Let us denote by W the region of the entire 
system that includes both the solute and solvent re- 
gions. Let us denote by V the region of solutes, or the 
cavity region, which is empty of solvent. We iden- 
tify the solute-solvent interface to be the boundary 
of the solute region V, and denote it by F = dV. We 
assume that the surface F consists of possibly many 
connected components, each of which is closed and 
continuous. 

For the cavity region V, we assign a volume exclu- 
sion function 



v{r) 



for f e V, 

1 else. 



Mathematically, this is the characteristic function of 
the solvent region W \ V which is the set of points 
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in W but not in V. The volume Vol [V] of the solute 
region V and the surface area Area [F] of the inter- 
face r can then be expressed as functionals of the 
volume exclusion function v{r) via 

Vol[V]= [ d^r^ [ [l-v{f)]d^r, 

Area[r]= [ dS ^ [ |Vw(f)|dV, 

where V = Vf is the usual gradient operator with 
respect to the position vector r and |Vw(f)| is the 
(5-function concentrated on the boundary T = dV 
of the cavity region V. The expression dS = 
Vw(f)|d^r can thus be identified as the infinitesi- 
mal surface element. We remark that, within the 
framework of the variational implicit solvent model 
[3 US I either the volume exclusion function v of 
the cavity or its boundary F can be used as the ulti- 
mate, direct variable of the solvation free energy of 
an underlying system. 

We assume that the position of each solute atom 
ri and the solute conformation are fixed. Thus, the 
solutes can be considered as an external potential 
to the solvent without any degrees of freedom. In 
this continuum solvent model, the solvent density 
distribution is simply = pov{'r), where po is the 
bulk density of the solvent. This means that we use 
a sharp interface approximation. 

B. Free energy functional 

For a given solvation system characterized by the 
cavity region V with its boundary F, the solute- 
solvent interface, the following ansatz of the Gibbs 
free energy was proposed in ^J, as a functional of 
the volume exclusion function v{ir) or its boundary 
F: 

G[v] = G™l[«] + GsurM + Gvdwb] + Gc\c[v] 

= P [l-v{f)]d^r+ j{f,r)\Wv{f)\d^i 
Jw Jw 

+ pa [ U{r}v{f)d^r + G,i,[v]. (1) 
Jw 

The first term, 

GvoiH = P f [1 - v{r)] dV = P Vol [V], (2) 
Jw 

proportional to the volume of V, is the energy of cre- 
ating a cavity in the solvent against the difference in 



bulk pressure P between the liquid and vapor phase, 
P = Pi ^ P^. In water at normal conditions or any 
fiuid close to phase coexistence, this pressure differ- 
ence is small and can often be neglected for solutes 
of microscopic size ('^nm). 
The second term 

GsurH= / 7(^r)|Vt;(r)|d3r= / l{?.T)dS 
Jw Jt 

(3) 

describes the energetic cost due to the solvent rear- 
rangement around the cavity, i.e., near the solute- 
solvent interface F, in terms of a function 7(f, F) 
with dimensions of free energy/surface area. This 
surface energy penalty is thought to be the main 
driving force behind hydrophobic phenomena 
It is a solvent specific quantity that also depends on 
the particular topology of the solute-solvent inter- 
face and varies locally in space p^ . 

The exact form of 7(r, F) is not known. The fol- 
lowing approximation based on a first-order curva- 
ture correction from scaled-particle theory ^2^ was 
made in 

7(r,F)=7o[l-2Ti/(f)], (4) 

where 70 is the constant solvent liquid-vapor surface 
tension for a planar interface, r is a positive constant 
often called the Tolman length [2l|, and 

H{r)^]^[n,{r) + K^(r)] 

is the local mean curvature in which ki (r) and K2 (r) 
are the local principal curvatures of the interface F. 

It has been shown in MD simulations that the sur- 
face tension 70 is the asymptotic value of the solva- 
tion free energy per unit surface area for hard spheri- 
cal cavities in water in the limit of large radii [10, 22]. 
In this system the Tolman length has been estimated 
to be of molecular size and has values of 0.7-0.9A. As 
its exact value is not known, the Tolman length may 
serve as the only fitting parameter in the variational 
continuum solvent model. The mean curvature H is 
defined only on the solute-solvent interface F. We 
have chosen the convention in which the curvatures 
are positive for convex surfaces (e.g., a spherical cav- 
ity) and negative for concave surfaces (e.g., a spher- 
ical droplet). 

We remark that, by the Hadwiger Theorem in dif- 
ferential geometry [23|, the geometrical part of the 
energy as a valuation of the closed surface F should 
have all the terms in Gvoi + Gsur (volume, surface 
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area, and surface integral of the mean curvature) 
plus an additional term of the surface integral of the 
Gaussian curvature 

K{r) = Ki{r)K2ir)- 

But, by the Gauss-Bonnet Theorem [23|, the Gaus- 
sian curvature is an intrinsic geometric property of 
the surface F, and its contribution to the free en- 
ergy is an additive constant. Therefore, it does not 
change our energy minimization process. We note 
that the Hadwiger Theorem was used in a general- 
ization of the classical theory of capillarity '25'| , and 
rece ntly in a morphometric approach to solvation 

ESUJ. 

The third term 



GvdwH = Po / Uir)v{f) &'r = po / C/(f)d-^r 
Jw JW\V 

(5) 

is the total energy of the non-electrostatic, van der 
Waals type, solute-solvent interaction given a sol- 
vent density distribution pQv{r). The potential 



N 



(6) 



is the sum of Ui that describes the interaction of 
the ith solute atom (with A'' total atoms) centered 
at fi with the surrounding solvent. Each term Ui 
includes the short-ranged repulsive exclusion and 
the long-ranged attractive dispersion interaction be- 
tween each solute atom i at position fi and a solvent 
molecule at f. Classical solvation studies typically 
represent the interaction Ui as an isotropic Lennard- 
Jones (LJ) potential. 



C/Lj(r) =4e 



(7) 
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(7) 



with an energy scale e, length scale cr, and center- 
to-center distance r. 

The last term Gcio[i'] is the electrostatic energy 
due to charges possibly carried by solute atoms and 
the ions in the solvent. In this work, we only consider 
nonpolar solutes. Therefore, we shall neg lect this 
term in what follows and refer to [1,0, mm and 
our forthcoming work [29l| for details. With this and 
considerations ([l])-([6]), we find for the final form of 



the nonpolar free energy functional 

G[v] = P Vol [V] + J 70 [1 - 2TH{f)] dS 

N 



Po 



J2 f U,{\f^n\)d'r, (8) 
■ , Jw\v 



where each interaction potential Ui {1 < i < N) has 
the form 



C. Free energy minimization 

Let Vmini'r) with its boundary r„iin be the ex- 
clusion function which minimizes the functional ([8]). 
Then, the resulting Gibbs free energy of the system 
is given by G'[winin]- The solvation free energy AG is 
the reversible work to solvate the solute and is given 

by 

AG = GKin] - Go, 

where Go is a constant reference energy which can 
refer to the pure solvent state and an unsolvated so- 
lute. The solvent-mediated potential of mean force 
along a given reaction coordinate x (e.g., the dis- 
tance between two solute centers of mass) is given, 
up to an additive constant, by w{x) — G[wmin], 
where Wmin('^) must be evaluated for every x. 

A necessary condition for T to be an energy mini- 
mizing solute-solvent interface is that the first varia- 
tion of the free energy functional ([5]) vanishes at the 
corresponding volume exclusion function v, i.e., 

SG[v 



Sv 



= 



(9) 



at every point of the boundary T. This energy vari- 
ation can be identified as a distribution over the in- 
terface r, and is given by p^[l3| 



SG[v] 
6v 



P + 2-fo[Hif)-rK{f)]-poUif). (10) 



The partial differential equation (PDE) determined 
by ([9]) and (fTO|) for the optimal exclusion function 
Wmin(^), or equivalently the optimal solute-solvent 
interface Fmin, is expressed in terms of pressure, cur- 
vatures, short-range repulsion, and dispersion, all 
of which have dimensions of energy density. It can 
be interpreted as a mechanical balance between the 
forces per surface area generated by each of the par- 
ticular contributions. A similar expression without 
the dispersion term was derived by Boruvka and 
Neumann within a generalization of classical cap- 
illarity (25.] . 
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III. THE LEVEL-SET METHOD 
A. Basics 

The starting point of the level-set method is to 
identify a surface T in three-dimensional space as 
the zero level-set of a function (f> — (f>{r) [l^, [l^, [l3| : 

r = {f: (j){f) = 0}. 

This means that the surface consists exactly of those 
points r at which the function (j) vanishes. This is 
in contrast with a parametric description of the sur- 
face r : r = f{a,(3) with a,P the parameters. The 
function = 0(r) is called a level-set function of 
the surface F. Clearly, the level-set function whose 
zero level-set represents the surface F is vastly non 
unique. 

The level-set function (p can be used to calculate 
many important geometrical quantities of the sur- 
face F. For instance, the unit normal vector n at the 
interface F, the mean curvature H, and the Gaus- 
sian curvature K can all be expressed in terms of 
the level-set function (h: 



1. 



H ^ -V -n, K ^n-adj {He{(t>))n, 

(11) 

where He{<j)) is the 3x3 Hessian matrix of the 
function whose entries are all the second order 
partial derivatives df^4) of the level-set function 0, 
and adj (-ffe((/))) is the adjoint matrix of the Hessian 
He{4>).^ 

Consider now a moving surface F — T{t) at time t. 
Let (f) = (f>{r, t) be a level-set function that represents 
the surface T{t) at time t. The basic idea is now 
to track the motion of the moving surface T{t) by 
evolving the level-set function (j)(r, t) and its zero 
level-set at each time t. The level-set function is 
determined by the so-called level-set equation, 



dt(f)-^v„\y(f)\ = 0, 



(12) 



where u„ is the normal velocity at the point f on the 
surface T{t). This normal velocity w„ = w„(r(i)) of 
each point r = r{t) on the surface F — T{t) at time 
t is defined by 



Vn = Vn{r{t)) 



dr{t) 
dt 



(13) 



The velocity is usually extended away from the sur- 
face so that the level-set equation (IT^ can be solved 
in a finite computational box. 



One of the major advantages of the level-set 
method is its easy handling of topological changes of 
surfaces during the surface evolution. For instance, 
the merge or break of bubbles can be captured by 
level-set calculations. This method is thus a perfect 
choice for capturing different kinds of stable solute- 
solvent interface structures. 



B. Normal velocity 

We apply the level-set method to evolve an initial 
interface to an equilibrium solute-solvent interface. 
This means that our level-set evolution is an opti- 
mization process rather than the real dynamics of 
the solvation system. We need to choose the gov- 
erning normal velocity of the interface in such a way 
that it will decrease the free energy during the sur- 
face evolution. As in common practice, we define 
the normal velocity of level-set to be the negative of 
the first variation of the Gibbs free energy: 



6G[v] 
Sv 



(14) 



By pU|) , the normal velocity is a function defined on 
F, and is given by 

vn = -P-2-/a[H{r)-TK{r)]+poUir}. (15) 

Here, we choose the unit normal n at F to point from 
the solute to the solvent region. 

Notice that the interface F — T{t), the volume 
exclusion function v = v(t), and the normal velocity 
Vn = Vn{t) all depend on the time t. This is not 
the time in the real dynamics but rather represents 
the state of optimization iteration. In particular, 
the normal velocity does not represent that of the 
interface evolution in real dynamics of the system, 
and hence can have non-physical units. It follows 
from the Chain Rule and the definition of the normal 
velocity (fT5)) that the time derivative of the Gibbs 
free energy G[w(t)] is 



dt 



GHt)] 



r(t) 



Sv 



dr{t) 





SG[v{t)] 


h(t) 


Sv 



dt 

2 



dS 



dS < 0. 



This shows that the normal velocity defined by 
decreases the energy. 
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With a given initial surface, we solve the level-set 
equation (jl2p each time step with the normal veloc- 
ity given by (fT5)) until a steady-state is reached. The 
steady-state solution gives a stable, energy minimiz- 
ing, solute-solvent interface. 

C. Implementation 

Our level-set algorithm consists of the following 
steps: 

Step 1. Input parameters and initialize an inter- 
face. The physical parameters include the pressure 
difference P, the macroscopic surface tension 70, the 
Tolman length t, the water density pq, the LJ energy 
and length parameters and Ui , and the coordinates 
of centers of all the fixed solute atoms i. An initial 
interface is defined by its level-set function; 

Step 2. Calculate the unit normal n, the mean 
curvature 7J, and the Gaussian curvature K by pT|) . 
Calculate the free energy using the formula ([5]); 

Step 3. Calculate the normal velocity Vn using the 
formula (jlSp . and extend the normal velocity w„ to 
the entire computational domain; 

Step 4- Solve the level-set equation p^ . As usual, 
this equation needs to be solved only locally near the 
interface F, since the value of (j) away from F will not 
affect the location of F; 

Step 5. Reinitialize the level-set function 0. The 
gradient of a solution to the level-set equation 
at certain time t can be sometimes too large or too 
small. This can lead to an inaccurate approximation 
of the normal n, the mean curvature -ff, and the 
Gaussian curvature K by (jlip . The reinitialization 
process uses the solution of the level-set equation 
(IT2t to obtain a new level-set function cj) that has 
the same zero level-set, i.e., the location F is not 
changed, and that keeps the gradient | V(/)| away from 
or from being too large; 

Step 6. Locate the interface F by the level-set 
function obtained in the previous step. Update the 
time step and go back to Step 2. 

There are two sources of instability in our level-set 
calculations. One is the Gaussian curvature term in 
the normal velocity p^ . Elementary calculations 
show that the motion by the combination of the 
mean and Gaussian curvatures, the H and K terms 
in (jl5[) . results in a parabolic equation of degener- 
ate type in certain parameter regimes. Specifically, 
one or two of the three eigenvalues of a matrix that 
defines the type of differential equations can become 
or even negative. In such a case, we add a small. 



positive constant to such eigenvalues so that the evo- 
lution is regularized. Such a regularization does not 
affect the final equilibrium solution. 

The other is the rapid change of values of the 
Lennard- Jones potential ([7]) when the distance is 
small. This can easily lead to large errors in the cal- 
culation of the free energy. To deal with this insta- 
bility, we have developed a two-grid algorithm. Our 
idea is to evolve the level-set function on a coarse 
grid and to calculate the energy on a fine grid. In 
doing so, we use interpolation and projection tech- 
niques to pass along information between the two 
grids. Numerical results show that this treatment 
works very well. 

IV. APPLICATIONS 

In this section, we report on our level-set calcula- 
tions of four nonpolar systems and compare the re- 
sults to available MD simulations using the SPC/E 
water model [s^. The four systems are two xenon 
atoms, two helical alkane chains, a single fuUerene 
Ceo, and two paraffin plates. 

In all of our calculations, we focus on water close 
to normal conditions (T — 298K and P=lbar) so 
that the pressure term in the free energy can be ne- 
glected. The other parameters are the macroscopic 
surface tension 70, the Tolman length r, and the wa- 
ter density po- The surface tension for SPC/E water 
has been calculated to be 70 = 72niN / m in agree- 
ment with the value of real water [31|. The den- 
sity is po = O.OSSA"^. The Tolman length has been 
estimated roughly for SPC/E to be r = 0.9A ^ 
and will serve as our only fitting parameter. All 
solutes are modeled by assemblies of identical and 
uncharged spheres interacting with the LJ potential 
([7]) with energy and length parameters e and a as 
summarized in Tab. I. These solutes are assumed 
to be in fixed configurations and have no degrees of 
freedom. 

For the helical alkanes and the Ceo, we perform 
additional MD simulations to provide data for the 
solvation free energy that is not available in litera- 
ture. The simulations are done using the MD simu- 
lation package DLPOLY [s^l in the NPT ensemble 
with up to iV = 800 SPC/E water molecules. The 
solvation free energy is calculated using standard 
thermodynamic integration ^33j procedures. Here, 
at least 15 simulation runs for different integration 
parameters A G [0, 1] per system are considered with 
lOOps equilibration time and Ins for gathering statis- 
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System 


e/fcsT 


a/A 


Xenon 


0.431 


3.57 


Helix (CH2) 


0.265 


3.54 


Ceo (C) 


0.158 


3.19 


Plate (CH2) 


0.265 


3.54 



TABLE L Investigated system LJ parameters: the atom- 
water LJ energy parameter e is in units of the thermal 
energy fcsT, and the atom-water LJ length a is in A. 



tics, where A corresponds to the scaled LJ length. 
The obtained ensemble-averages were interpolated 
by an Akima spline and the resulting curve inte- 
grated to get the solvation free energy. 

In the level-set method, we usually start with a 
large spherical surface in a cubic box that encloses 
all the fixed solute atoms and then evolve the sur- 
face to minimize the nonpolar solvation energy. The 
box length is between 15-25A and a grid size of 
50'^ or 100'^ bins is used depending on the solute 
size and desired computational speed and accuracy. 
Finite-size corrections to the dispersion part of the 
energy are considered by integrating the long-ranged 
LJ interactions over a homogeneous water distribu- 
tion beyond the box up to infinity. The computa- 
tional time of a single level-set minimization takes 
several hours on a single-processor computer but sig- 
nificantly depends on the choice of the initial con- 
figuration. Independent of the latter, we find that 
the system always converges to a stable configura- 
tion corresponding to an energy minimum. We have 
not investigated whether the systems exhibit multi- 
ple energy minima and postpone this rather complex 
issue to a future study. 



A. Two xenon atoms 

In this example, we investigate the performance of 
our method for a simple system of microscopic non- 
polar solutes that consists of two xenon (Xe) atoms. 
For that, we fix the Xe atoms in a center-to-center 
distance d and calculate the solvent-mediated inter- 
action wld) between them, which is basically the 
solvation energy in dependence of the coordinate d. 
The total potential of mean force (pmf) W{d) is 
the sum of the solvent-mediated part and the intrin- 
sic Xe-Xe vacuum interaction which can be modeled 
also by an LJ interaction. The water-xenon LJ in- 
teractions we use for our level-set calculations are 



taken from Paschek [3J| who accurately calculated 
the pmf in MD simulations. For our comparison, 
we find that Tolman lengths between r = 0.9 and 
l.OA give reasonable agreement when compared to 
the MD results, close to the value of 0.9A estimated 
by previous MD simulations. 



1 1 1 1 1 1 1 , 1 , 1 1 1 1 1 1 1 
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FIG. 1: Comparison of the interaction between two 
xenon atoms from level-set and MD calculations. The 
solid line and dotted line with symbols (circles) are the 
solvent-mediated interaction w{d) by the MD and level- 
set calculations, respectively. The inset displays the full 
pmf W(d) where the vacuum Xe-Xe LJ interaction is 
added to the solvent-mediated part. 

In Fig. 1, we plot the solvent-mediated part wld) 
and the total pmf W{d) of our level-set calcula- 
tion compared to the MD simulation reported by 
Paschek [3J| for a Tolman length r = 0.95A. Also 
plotted are interface images from the level-set so- 
lution at three selected distances. For small sep- 
arations d < 4.5A, the solution gives two over- 
lapping spheres and the solvent-mediated interac- 
tion is attractive in agreement with the MD simula- 
tions. The attraction comes from a smaller water- 
accessible surface due to the overlap of spheres and 
the gain of resulting interfacial energy. At sepa- 
rations d ~ a maximum of about Ifc^T in 
height occurs in both MD and level-set calculations, 
also in good agreement with each other. In the con- 
tinuum approach, this desolvation barrier is implic- 
itly accounted for by the unfavorable LJ interaction 
when replacing water molecules adjacent to the first 
Xe atom by the second Xe atom. For separations 
d > 6A, the interface breaks (water penetrates) and 
the stable level-set solution corresponds to two iso- 
lated spheres. The shallow oscillations in the MD 
curve for d > 7.5A are due to explicit water struc- 
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turing around the Xe atoms and are not captured 
by our continuum method. OveraU, however, we can 
judge that the agreement is good and the dominant 
features of the pmf are well described by our macro- 
scopic method using just one fit-parameter which is 
close to previous estimates. 



B. Helical alkane chains 

In order to check how our level-set method per- 
forms for larger and more complex shaped molecules, 
we study model helical alkane chains that are as- 
sembled by CH2 atoms modeled by the OPLS-AA 
force field [IHl; see Tab. I. We investigate the solva- 
tion of two different configurations, one more loosely 
packed involving 32 atoms (alkane A) and the other 
one more tightly packed using 22 atoms with hardly 
room for water in the helical core (alkane B). 

The results are plotted in Fig. 2 which include a 
comparison of our level-set calculation and a typical 
SASA type surface constructed by taking the enve- 
lope of all the LJ spheres. Though they look quite 
similar, our level-set result leads to a much smoother 
solute-solvent interface, a result from the minimiza- 
tion of interface area based on the energy functional. 




FIG. 2: A comparison of the level-set (left) and SASA 
(right) description of helical polymer chain A with 32 
atoms (above) and B with 22 atoms (below), respec- 
tively. 



From our MD simulations, we estimate solvation 
free energies of AG ~ Gfe^T and AG ~ IksT for 
the alkanes A and B, respectively. The energies are 
positive and small and have the same order of magni- 
tude as the solvation energy of small alkanes f36l.l37t. 
These relatively small free energies are a well-known 
consequence of enthalpy-entropy compensation in 
the solvation of nonpolar solutes |36l.[37|. It has been 
shown that AG can be decomposed in solute-solvent 
enthalpy A[/uv and solute-solvent entropy TASuv, 
while the solvent-solvent enthalpy and entropy ex- 
actly cancel each other [H, H^] ■ For nonpolar solutes 
the enthalpic part is the (ensemble-averaged) total 
solute-solvent LJ interaction which we estimate by 
our MD simulations to be a large AUuv — -lOTfcsT 
and -b'aksT for the alkanes A and B, respectively, 
indicating solute-solvent entropies of the same order 
of magnitude. 

In the level-set method, we can reproduce the to- 
tal solvation free energy for both alkanes within 5% 
by using a Tolman length of r = 1.3 A. In our implicit 
model, we can identify the solute-solvent enthalpic 
part A[/uv by the LJ term Gvdw =: Gvdwbmin] 
and the entropic part —TASuv by the interfacial 
term Gsur —■ Gsui[^'min] in the free energy func- 
tional; wc find large values of Gvdw = — OSfc^T and 
Gsur = lOlksT for helix A, and Gvdw = -^SksT 
and Gsur — ^bksT for helix B, close to the values 
obtained by the MD simulations. We note that small 
variations of the Tolman length r have a negligible 
effect on the enthalpic part of the free energy while 
its influence on the entropic (interfacial) part is no- 
ticeable. For example, we find Gsur — 63fcsT for he- 
lix B when using r = 1.1 A instead of Gsur — 55k bT 
for r = 1.3 A. As the total solvation free energy is 
a difference of two large quantities, the interfacial 
part of the free energy functional ([1]), in particu- 
lar the curvature correction, has to be reconsidered 
carefully. 

In Fig. 3, we show the same two helical chains but 
using a color code that indicates the value of mean 
curvature at each point of the solute-solvent inter- 
face. The curvature varies between positive and neg- 
ative values, showing as well concave parts (blue). 
The highest positive curvature (deep red) is roughly 
given by the inverse of the length a of one LJ sphere. 
Note that the concave parts of the loosely packed 
helix are within the helical core in contrast to the 
convex outer parts. This qualitatively different hy- 
dration of the helix depending on the local (convex 
or concave) geometry is in line with structural stud- 
ies of water at protein- water interfaces . 
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FIG. 3: Level-set calculations of helical molecules. Color 
represents mean curvature. The units of mean curvature 
are in A~^. 



These examples show that complex shapes with 
varying interface curvature, as typical in protein 
structures, can be efficiently tackled by the vari- 
ational implicit solvent model in conjunction with 
the level-set method. The solute-solvent interface 
and the resulting energies seem to be well described 
by our methods, in particular when regarding the 
fact that the solvation free energy is a difference of 
large entropic and enthalpic contributions which of- 
ten leads to a large error in the calculation of the 
total solvation free energy. 



C. A single molecule of fullerenes Ceo 

Another interesting molecule which can be rea- 
sonably modeled as a nonpolar entity is the Cgo 
molecule, the buckyball. The C-atom-water LJ pa- 
rameters are taken from [s^ and are also shown in 
Tab. I. The experimental solvation free energy in wa- 
ter is typically close to zero 40] in agreement with 
our MD simulations which yield AG ~ — Ifc^T. Pre- 
vious more empirical implicit solvent studies show 
that the large interfacial energy penalty is more than 
compensated by the strong dispersion attraction be- 
tween the water and the tightly packed carbon shell 
resultin g iii a small and negative total solvation free 
energy [40|. A related unusual repulsive solvent- 
mediated interaction between two Cm molecules has 
been observed in MD simulations [41[. 

Our level-set result for the equilibrium interface is 
shown in Fig. 4, featuring a smooth soccer-like sur- 
face. Using a Tolman length of r = I.SA, we can 
reproduce the MD solvation energy within 5%. A 
separate analysis of the interfacial energy part and 
the dispersion part shows that a large AGvdw ~ 
—AQksT is gained from dispersion attraction while 
AGsur = 'iSksT surface energy is paid. Our MD 



simulations confirm this cancellation of energy con- 
tributions quantitatively by showing enthalpic and 
entropic contributions of ~ — SOfcsT and ~ 49kBT, 
respectively. 




FIG. 4: A comparison of the level-set (left) and a SASA 
type (right) description of a Ceo buckyball. 

In Fig. 5, we show the same Ggo molecule obtained 
by our level-set calculation with a color code that in- 
dicates the value of mean curvature of the interface. 
In contrast to the helical molecules, we find only 
convex curvatures varying from zero to roughly the 
inverse of the LJ size of one C atom. The curva- 
ture distribution displays the typical five and sixfold 
structure of the Cgo molecule. 




FIG. 5: The level-set calculation of the Ceo molecule 
with color indicating the value of mean curvature. The 
units of mean curvature are in A~^. 

As in the previous case of the two helical alka- 
nes, this example demonstrates that the subtle bal- 
ance between interface (entropic) and dispersion (en- 
thalpic) terms, and thus the correct interface loca- 
tion and curvature are crucial for an accurate de- 
scription of solvation free energies, and are well cap- 
tured by our methods. 
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D. Two parallel nanometer-sized paraffin 
plates 

In our last example, we consider the strongly hy- 
drophobic system of two parallel paraffin plates as 
investigated in the explicit water MD simulations 
by Koishi et al. Each plate consists of 6 x 6 

fixed CH2-atoms with atom-water LJ parameters 
from the OPLS-AA force field, see Tab. 1, and has 
a square length of ~3nm. The two plates are placed 
in a center-to-center distance d and different separa- 
tions are investigated. Koishi et al. observed a clear 
dewetting transition (vapor bubble or "nanobubble" 
formation) for distances d < 15 A accompanied by a 
strong attractive interaction energy of the order of 
tens of ksT. The pmf is shown in Fig. 6 together 
with the solution of our variational implicit solvent 
model and level-set snapshots of the equilibrium in- 
terface at selected distances. 




FIG. 6: Comparison of the pmf from level-set (circles) 
and MD calculations (dashed line) for the two-plate sys- 
tem. Three level-set interfaces of the system are shown 
at their respective separations. 

As can be seen in Fig. 6, we have almost per- 
fect agreement with the MD simulation results. We 
find that the curve hardly depends on the particular 
choice of the Tolman length (t — 0.9A here). The 
reason is that the average radius of mean curvature 
for this large-scale example is much bigger than the 
typical value of the Tolman length and the theory 
basically becomes fit-parameter free. As another im- 
portant result, our variational implicit solvent model 
captures the nanobubble formation, see the interface 
snapshots in Fig. 6 for c? < IsA, a consequence of 
the energetic desire of the system to minimize the 



total interface area. For larger distances the inter- 
face breaks, i.e., the bubble ruptures and the equilib- 
rium interface is given by two isolated plates having 
almost no mutual interaction. 

We note here that for distances d < 15 A where 
bubble formation takes place, the free energy of the 
system might have a second, local minimum corre- 
sponding to the wetted case with two isolated plates 
[10]. The existence of two local minima shows the 
possibility of hysteresis in bubble formation as ob- 
served in the MD simulations [42]. 

V. CONCLUSIONS AND OUTLOOK 

We have developed a level-set method for numer- 
ically capturing stable solute-solvent interfaces for 
complex shaped nonpolar molecules in three dimen- 
sions based on a recently developed variational im- 
plicit solvent model. This method evolves an ar- 
bitrary initial surface to decrease the solvation free 
energy until a stable equilibrium solute-solvent in- 
terface is reached. 

In implementing our numerical method, we have 
developed a regularizing technique to stabilize the 
surface evolution governed geometrically by a com- 
bination of both the mean curvature and Gaussian 
curvature. We have additionally developed a stable, 
two-grid method for the calculation of the total free 
energy. 

We have applied our method to several nonpolar 
systems. Our numerical results show good agree- 
ment with MD simulations with reasonable interface 
shapes and curvatures. In particular, due to the 
numerical robustness and ability of handling topo- 
logical changes, our method captures the volume 
fusion/break and nanobubble formation/rupture in 
the two xenon system and the strongly hydrophobic 
two-plate system, respectively. 

A comment is needed for the treatment of large 
concave curvatures which can occur locally in micro- 
scopic systems. In the two Xe atom example imme- 
diately before break-up, a singularity near the neck 
of the two spheres develops and can artificially in- 
crease the energy. We find that the energy is lowered 
if we renormalize the Tolman length for very large 
mean curvatures. Mathematically, it is known that 
even the motion of surface solely by mean curvature 
can induce the neck-pinch singularity that we see in 
our level-set calculations ^3|. Therefore, it remains 
to be further investigated how the singularity should 
be treated in the free energy definition and level-set 
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calculation. 

We emphasize that our level-set implicit solvent 
calculations are one or two orders of magnitude 
faster than explicit MD simulations. Our method 
does require the solution of the level-set equation 
([T^ with the normal velocity in each time step. 
Compared with a SASA type approach, this is an 
"extra" work, but can be done efficiently. For in- 
stance, if we start with an initial surface that is close 
to an equilibrium one, our calculations can be much 
more efficient. With a good initial guess and a rea- 
sonable grid size, we need about 15-20 minutes to 
run our code for the calculation of a helical polymer 
chain. When the level-set method is applied to real 
dynamics calculations using continuum models, the 
relaxation of the surface to a complete equilibrium is 
not necessary, and only a few iterations are enough 
to update an interface. Therefore, in such a case, the 
level-set calculation can be compatible with a SASA 
type method in terms of efhciency. 

Our level-set method can be used for the force 
calculation of an underlying solvation system. In 
principle, forces are obtained by the gradient of the 
free energy with respect to some spatial coordinates. 
One such coordinate is the geometrical location of an 
equilibrium solute-solvent interface. Our artificial or 
optimization normal velocity is exactly the effective 
force with respect to such a coordinate. This re- 
places the calculation of surface area in a SASA type 
method. Within the framework of level-set method, 
we can also evaluate forces between the solute atoms 
that can be used as input to Brownian dynamics 
computer simulations. This is an important issue 
that we are still investigating in details. 

Finally, let us comment on further developments 
which will be crucial for a complete implicit sol- 
vation description of large biomolecules by our 



method. We are currently developing a level-set 
method for capturing numerically stable solute- 
solvent interfaces using the Gibbs free energy ^ 
that includes the electrostatic contribution of the 
polar groups of the solutes. This method couples 
the presented level-set method to a finite-difference 
based solver for the Poisson-Boltzmann (PB) equa- 
tion, a typical approach for the implicit treatment 
of electrostatics in solvated molecular systems |6l, ^ . 
In [29] , we derive the variation of the free energy ([1]) 
including electrostatics on the PB level with respect 
to the change of the interface P. This will be used 
to define the normal velocity f„ similar to that in 
(fT4)) but extended to local electrostatic pressures. 
Currently we are also developing fast and optimized 
level-set methods for the handling of large systems 
that can have a few thousand atoms in solutes and 
that can allow solute atoms to freely move around 
in the optimization process. Another challenge is 
to extend the level-set treatment to predict the 
correct dynamics of evolution of the interface [4^ 
and account for interface fluctuations. 
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